home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_12_01 / letters / ratint.c < prev    next >
Text File  |  1993-06-07  |  1KB  |  42 lines
  1.                          Listing 1
  2.              Function RATINT} for the RPFT Code
  3.  
  4. /******************************************************
  5.  * RATINT - Diagonal rational function interpolation in
  6.  * the arrays xa[1..n] and ya[1..n].
  7.  *****************************************************/
  8. void ratint(double xa[], double ya[], double *c,
  9.         double *d, int n, double x, double *y)
  10. { int m,i,ns=1;
  11.   double w,t,hh,h,dd;
  12.   static double miny=1.e99;
  13.  
  14.   if (miny>1.e90) for (i=1; i<=n; ++i)
  15.                          if (ya[i]<miny) miny=ya[i];
  16.   hh=fabs(x-xa[1]);
  17.   for (i=1;i<=n;i++)
  18.   { h=fabs(x-xa[i]);
  19.     if (h == 0.0)  { *y=ya[i]; return; }
  20.     else  if (h < hh) { ns=i; hh=h; }
  21.     c[i]=ya[i]-miny; d[i]=ya[i]-miny+1.e-50;
  22.   }
  23.   *y=ya[ns--] - miny;
  24.   for (m=1;m<n;m++)
  25.   { for (i=1;i<=n-m;i++)
  26.     { w=c[i+1]-d[i] ; h=xa[i+m]-x; t=(xa[i]-x)*d[i]/h;
  27.       dd=t-c[i+1];
  28.       if (fabs(t)>1.e15) 
  29.         fprintf(stderr,"Probable loss of accuracy in "
  30.           "RATINT. fabs(t) > 1.e15 for X = %.8G\n",x); 
  31.       if (dd == 0.0)
  32.       { fprintf(stderr,"Error in routine ratint. The "
  33.            "function may have a pole at x=%.8G\n",x);
  34.         exit(1);
  35.       }
  36.       dd=w/dd; d[i]=c[i+1]*dd; c[i]=t*dd;
  37.     }
  38.     *y += (2*ns < (n-m) ? c[ns+1] : d[ns--]);
  39.   }
  40.   *y += miny;  return;
  41. }
  42.